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Abstract 

We describe a Fourier accelerated Hybrid Monte Carlo algorithm 
suitable for dynamical fermion simulations of non-gauge models. We 
test the algorithm in supersymmetric quantum mechanics viewed as 
a one-dimensional Euclidean lattice field theory. We find dramatic 
reductions in the autocorrelation time of the algorithm in comparison 
to standard HMC. 
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Introduction 


Dynamical fermion algorithms play a crucial role in the simulation of lattice 
field theories. The favorite algorithm for an even number of fermion species 
is the Hybrid Monte Carlo HMC algorithm |j|]. In this paper we introduce 
and test an improved version of this algorithm in the case of supersymmetric 
quantum mechanics. We show that the improved algorithm is far superior to 
the usual HMC procedure in combating the effects of critical slowing down. 

We first discuss the usual HMC algorithm and show how it may be gen¬ 
eralized to allow for Fourier acceleration. The idea is closely related to the 
usual Fourier acceleration used for Langevin simulations - both the fields and 
their conjugate momenta are evolved in momentum space throughout an in¬ 
dividual HMC trajectory - the crucial improvement being to choose a wave¬ 
length dependent timestep. The momentum dependence of this timestep 
is chosen so as to render the magnitude of the resulting update of a given 
Fourier component insensitive to its wavevector index. We test this algo¬ 
rithm in supersymmetric quantum mechanics treated as a Euclidean lattice 
theory. This model is formulated on the lattice in such way as to leave intact 
an exact subgroup of the continuum supersymmetry. We are able to demon¬ 
strate that the improved algorithm reduces the dynamical critical exponent 
to values close to zero on lattices as large as L = 256 and correlation lengths 


^ ~ 16. 


Hybrid Monte Carlo algorithm. 


We will be concerned with simulations of models involving scalar and fermion 
fields. The typical action we will be discussing takes the general form: 


S{x, Ip, Ip) = Sb{x) + 'ipM{x)ip 


( 1 ) 


containing a bosonic field x and a Dirac field ip defined on a lattice in 
Euclidean space. The fermion matrix M will contain lattice derivative terms 
together with couplings to the bosonic field x. The partition function for 
this system is then just 



( 2 ) 


In order to simulate this action we first replace the fermion field by a 
bosonic pseudofermion field (p whose action is 
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This is an exact representation of the original boson effective action in the 
case where det M > 0 (which will be the case for SUSY QM). The resultant 
nonlocal action S (x, (j)) can be simulated using the Hybrid Monte Carlo 
(HMC) algorithm |jl]. In the HMC scheme momentum fields {p, vr) conjugate 
to (x, 4>) are added and a Hamiltonian H constructed from the original action 
and additional terms depending on the momenta: 

H = S + AS 

= {Pi + ^i) 

i 

The corresponding partition function Z' = e~^ is, up to a constant factor, 
identical to the original Z. The augmented system {x,p,4>,tt) is now as¬ 
sociated with a classical dynamics depending on an auxiliary time variable 


dx 

dp 

dH 

'm ’ 

m ^ 

dx 

d<j) 
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dH 


'dt ^ 

dcj) 


Introducing a finite time step At allows us to simulate this classical 
evolution and to produce a sequence of configurations {x{t),(p{t)). 

If At = 0 then H would be conserved along such a trajectory. As At 
is finite H is not exactly conserved. However a finite length of such an ap¬ 
proximate trajectory can be used as a global move on the fields (x, (j)) which 
may then be subject to a metropolis step based on AH. Provided the dis¬ 
crete, classical dynamics is reversible and care is taken to ensure ergodicity 
the resulting move satisfies detailed balance and hence this dynamics will 
provide a simulation of Z' and hence also of the original partition function 
Z. 

Ergodicity is taken care of by drawing new momenta from a Gaussian 
distribution after each trajectory. The reversibility criterion can be satisfied 
by using a leapfrog integration scheme. Its general form for a field x with 
associated momentum p can be written as 

(St)'^ 'T 

Xr{St) = Xr(0) -|-(ItArr'Pr'(O) “I-^— Arr" A^Uy, Fr'{Q) 

Pr{6t) = p,{0)+^-^Al,{Fr>{0)+Fr>{dt)) (3) 

where F = —dH/dx is an associated force and A is an arbitrary matrix. 
Notice our notation - the fields are indexed by a integer vector giving their 
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lattice position. This update indeed satisfies the time reversibility criterion: 
inverting the sign of momentum as Pni^t) = —pniSt) and updating the 
resultant configuration leads to = 4>n{0), Pn(25t) = —pn{0). 

To implement our improved algorithm we choose a matrix Ay-j./ whose 
elements depend only on the difference of the lattice vectors between two 
sites r and r'. Ay.y.i = o (r — r') 

We then take the lattice Fourier transform of equation 3. to arrive at 
an update equation for the Fourier components and py. The Fourier 
components of the field Xy are given by the discrete transform 

= (4) 

with similar expressions for pk and Ok- Because of the lattice convolution 
theorem the dependence of this momentum space update on the Fourier 
coefficients Ok may then be completely absorbed into the definition of a 
wavelength dependent time step Jtk = The final update equations 

take the form 

Xk(t + 5tk) = Xk(t) +^tkPfc(t) + Fk(f) 

Pk{t + Stk) = Pk(i) + ^ (Tk(t) + Tk(t + 5tk)) (5) 

The final step in applying Fourier acceleration to this update is to choose 
the acceleration kernel Ok so as eliminate the wavelength dependence of 
the update in the free theory [p. For the update of a bosonic field this 
implies that we should utilize the square root of the momentum space lattice 
propagator. Such a propagator contains a mass parameter ruacc which can 
be tuned to optimize the update. In the free theory it should clearly be set 
to the bare lattice mass, but in a general interacting theory it is left as a 
parameter. In this paper we provide evidence that the autocorrelation time 
is optimized by setting m-acc to the approximate position of the massgap. 
The fermionic kernel is then chosen to be the inverse of the bosonic kernel. 

Notice the occurrence of the square root of the momentum space lattice 
propagator here - rather than the propagator itself which would be usual in 
Fourier accelerated Langevin algorithms - this reflects the fact that here the 
HMC update corresponds to a discrete second order differential equation in 
time unlike the first order Langevin equation. This choice of the square root 
propagator was not made in |^] and may explain, in part, why acceleration 
led to only small reductions in the autocorrelation time. 
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Model 


A simple way to demonstrate the effectiveness of this update is to consider 
snpersymmetric quantum mechanics Q| which contains a real scalar field 
X and two independent real fermionic fields 'll: and defined on a one¬ 
dimensional lattice of L sites with periodic boundary conditions imposed on 
both scalar and fermion fields: 

1 1 

5 " = — ^ {DijXj + Pi) {DijXj + Pi) + — ^ ipi (^Dij Pj^j^ ipj ( 6 ) 
ij ij 

For our simulations the quantity Pi and its derivative are chosen as: 

Pi = J2 

3 

Plj = Kij + 3gXi6ij 

The symmetric difference operator Dij and the Wilson mass matrix Kij are 
defined as: ^ 

Pij — 2 

Kij — nidij — “ISij) 

Here dimensionless lattice nnits are m = ruphysa, g = 5physa^ and x = 
a~2Xphys and the discrete momenta take on the values 2'Kk/L with k = 
0...L — 1. Notice that this model employs a non-standard boson action 
containing not the usual scalar lattice Laplacian but the square of the sym¬ 
metric difference operator. This is done in order to treat the fermions and 
bosons in a symmetric manner - indeed because of this the action is invariant 
under a single SUSY-like symmetry. 

5xi = 'il^ii 

= {DijXj + Pi)C 
Si/ji = 0 

Doubles in both bosonic and fermionic sectors are eliminated by means of 
the Wilson mass term K. The physics results from this study were published 

in ||. 

For bosonic and pseudo-fermionic field updates respectively we use the 
following timesteps which are simple inverses of each other. 

= At{macc + 2)/^sin^(27r/c/L) -|- (mace + 2sin^(7rA:/L))2 
= A.t\Jsw?{2TTk/L) -|- {niacc + 2 sin^(7rfe/L))2/(mace + 2) 
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Autocorrelation time 


Suppose that Q is some observable and Qt is a measurement of Q in con¬ 
figuration corresponding to Monte-Carlo time t. Then the autocorrelation 
function is defined as: 

< QoQt > - < Q >2 

Typically the function c{t) is approximately exponential c{t) = where 
r is the autocorrelation time - the time between decorrelated configurations. 
Clearly we can define r also from the relation 


r = 



e = ^2, c{t) 

t 


and this yields a robust way to estimate r even when c{t) is not exactly 
an exponential - measured this way it is sometimes referred to as the inte¬ 
grated autocorrelation time. It also provides a convenient way to find the 
autocorrelation time from a sequence of N Monte Carlo measurements 


^ ^ ^ < Q2 > _ < Q >2 

where we measure the autocorrelation function only out to M steps (we 
must choose M >> r). 

The efficiency of different simulation algorithms can be described in 
terms of the dynamic critical exponent z: 


r ~ (1/Mi,ttice)" = (i/Mphysical)" ~ 

An algorithm which resulted in z = 0 is said to suffer no critical slowing 
down while most local algorithms such as conventional HMC typically yield 
z ~ 2. 

Fig. I represents the autocorrelation time r as a function of mace for 
lattice SUSY QM with m = 10, g = 100 for a lattice with L = 64 sites. No¬ 
tice that the dimensionless parameter characterizing the interaction strength 
g/w? is unity in this case - we are in a strong coupling regime. The massgap 
of the theory for these values of the bare parameters can be calculated via 
Hamiltonian methods and yields a value rUg ~ 16.0. 

For large mace ^ oo the update can be shown to reduce to the usual 
local HMC algorithm and, as the plot makes clear, leads to a maximal value 
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Figure 1; Inr vs In ruacc for Supersymmetric Quantum Mechanics with 
L = 64, m = 10, <7 = 100. 


for T. Indeed the minimal value of r is can be seen to be achieved by setting 
TTT-acc approximately equal to the massgap of the model In mace ~ 2.7. In the 
vicinity of this point, one can see that the minimum autocorrelation Tmin is 
more than an order of magnitude smaller than the r achieved by the usual 
local algorithm for this lattice size. 

To derive a dynamic critical exponent we have also compared the ded- 
pendence of the autocorrelation time on lattice size L both for the usual 
local HMC and the Fourier accelerated algorithm with macc = 15.0. For 
the case macc = oo fig shows a plot of In r against In L together with a 
straight line fit yielding 2 ; ~ 1.7. The plot makes it clear that this is a lower 
bound for 2 ; so that in all likelihood 2 ; achieves a value of 2 ; = 2 for very large 
lattices. This is to be contrasted to fig which shows the same plot for 
T^a.cc = 15. A fit in this case yields a value of 2 ; consistent with zero. For the 
largest lattice size L = 256 we see that more than two orders of magnitude 
decrease in r are gained by use of the accelerated algorithm. 
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Figure 2: In r vs In L for Supersymmetric Quantum Mechanics with mace ~ 
oo, m = 10, = 100. Straight line fit yields z = 1.7 

Discussion 

We have shown that the use of Fourier acceleration can yield large payoffs 
in the simulation of non-gauge field theories with dynamical fermions. We 
have presented results which support this conclusion for supersymmetric 
quantum mechanics in which the HMC algorithms can be pushed to large 
[L = 256) lattice size and correlation length ~ 16). We have obtained 
similar, though less quantitative, results for two-dimensional Wess-Zumino 
models |^. 
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Figure 3: In r vs In L for Supersymmetric Quantum Mechanics with mace = 
15, m=10, g = 100. 
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